Data

Para determinar las imagenes de la zona afectada usamos los vectores provisto por la plataforma geojson.io.

Las imagenes satelitales las descargamos desde la plataforma de sentinel, via la api sentinelsat.

Las imágenes corresponden a los dias

  • 20200815: _S2A_MSIL2A_20200815T141741_N0214_R010_T20JLM20200815T182929.zip
  • 20200817: _S2B_MSIL2A_20200817T141049_N0214_R110_T20JLM20200817T163953.zip
  • 20200820: _S2B_MSIL2A_20200820T141739_N0214_R010_T20JLM20200820T182855.zip
  • 20200822: _S2A_MSIL2A_20200822T141051_N0214_R110_T20JLM20200822T182240.zip
# from http://geojson.io/#map=11/-30.5995/-64.2055
vector_file='data/villa_albertina.geojson'
# images from https://scihub.copernicus.eu using sentinelsat api
zip_20200815='data/S2A_MSIL2A_20200815T141741_N0214_R010_T20JLM_20200815T182929.zip'
zip_20200817='data/S2B_MSIL2A_20200817T141049_N0214_R110_T20JLM_20200817T163953.zip'
#
zip_20200820='data/S2B_MSIL2A_20200820T141739_N0214_R010_T20JLM_20200820T182855.zip'
zip_20200822='data/S2A_MSIL2A_20200822T141051_N0214_R110_T20JLM_20200822T182240.zip'

Productos y temporales

#collapse
class bands_etl(object):
    """
        
    """
    
    @staticmethod
    def get_ndi(band_1,band_2,weight=[1.0,1.0],dtype=rio.float32,nodata=-10000):
        ndi = (weight[0]*band_1.astype(float)-weight[1]*band_2.astype(float))/(band_1*weight[0]*+band_2*weight[1])
        ndi_val=np.where((band_1*weight[0]*+band_2*weight[1])>0,ndi,nodata)
        
        return ndi_val.astype(dtype)
    
    @staticmethod
    def get_band_paths(path_,pattern_=['*_10m.jp2']):
        """
        """
        names_bands={}
        all_names_bands={}
        for root, _, _ in os.walk(path_):
            for p in pattern_:
                ret=glob.glob(root+os.path.sep+p)
                if len(ret)>0:
                    j=0
                    for r in ret:
                        full_path=os.path.join(root, r)
                        name=os.path.basename(full_path)
                        band=name.split('_')[-2]
                        c={band:full_path}
                        names_bands.update(c)
                        all_names_bands.update({p+'_'+str(j):full_path})
                        j+=1
                else:
                    pass
        return names_bands,all_names_bands
    
    @staticmethod
    def basic_resampling(input_raster,output_raster,resampling_method='cubic',cell_size=20):
        """
        :param input_raster: file to be resampled
        :param output_raster: file to be written
        """
        warp_opts = [
            "-r",
            resampling_method,
            "-tr",
            str(cell_size),
            str(cell_size),
        ]
        #
        gdal.UseExceptions()
        ds = gdal.Warp(output_raster, input_raster, options=warp_opts)  # noqa
        del ds
    
    def set_band_names(self,band_dict):
        """
        :param band_dict: dict with band names and paths
        """
        self.band_names_=band_dict
        
    def get_band(self,b,open_flag):
        """
        """
        if open_flag:
            band_=rio.open(self.band_names_[b])
        else:
            band_=b
        return band_
    
    def set_mask_from_gpd(self,mask_gpd_geometry):
        """
        """
        self.mask_gpd_geom_=mask_gpd_geometry
    
      
    def get_clipping_params(self,band_init,mask_opts={'crop':True},open_flag=True):
        """
        """
        band_=self.get_band(band_init,open_flag)
        out_image_init, out_transform_init = mask(band_, self.mask_gpd_geom_,**mask_opts)
        self.image_shape_0_=out_image_init[0].shape[0]
        self.image_shape_1_=out_image_init[0].shape[1]
        self.transform_=out_transform_init
        self.dtype_=out_image_init[0].dtype
        self.crs_=band_.crs
        
    def bands_clipped_to_tiff(self,bands_list,tiff_name_='bands.tiff',driver_='Gtiff',open_flag=True,scale_band=False):
        # Create an RGB image 
        if scale_band:
            dtype_=rio.float32
        else:
            dtype_=self.dtype_
        
        
        with rio.open(tiff_name_,'w',driver=driver_, width=self.image_shape_1_, height=self.image_shape_0_, 
              count=len(bands_list),crs=self.crs_,transform=self.transform_, dtype=dtype_) as img:
    
            for i,b in enumerate(bands_list,start=1):
                band_=self.get_band(b,open_flag)
                out_image_init, out_transform_init = mask(band_, self.mask_gpd_geom_,crop=True)
                if scale_band:
                    res=out_image_init[0]/out_image_init[0].max()
                    img.write(res.astype(dtype_),i)
                else:
                    img.write(out_image_init[0],i)
                #
                band_.close()
                out_image_init=None
                out_transform_init=None
            img.close()
    
    
    def ndi_to_tiff(self,band1_,band2_,weight=[1.0,1.0],tiff_name_='ndi.tiff',driver_='Gtiff',open_flag=True,nodata=-10000):
        # Create an NDVI image 
        with rio.open(tiff_name_,'w',driver=driver_, width=self.image_shape_0_, height=self.image_shape_1_, 
              count=1,crs=self.crs_,transform=self.transform_, dtype=rio.float32,nodata=nodata) as img:
            band2_init=self.get_band(band2_,open_flag)
            band1_init=self.get_band(band1_,open_flag)
            band2, _ = mask(band2_init, self.mask_gpd_geom_,crop=True)
            band1, _ = mask(band1_init,  self.mask_gpd_geom_,crop=True)
            band2_init.close()
            band1_init.close()
            img.write(self.get_ndi(band1,band2,weight,nodata=nodata))
            img.close()
            
            
def write_array_single(
    input_filename,
    output_filename,
    array,
    array_dtype=rio.uint8,
    nodata_=0,
    attribs_source_="profile",
):
    """
    :param input_filename: input filename (vrt,tiff,etc) to be taken as src for profile
    :param output_filename: output filename
    :param array: numpy array to be written
    :param array_dtype: rasterio or numpy dtype
    :param nodata_: nodata to be written as integer (e.g.:-10000)
    :param attribs_source_: meta or profile
    """

    with rio.open(input_filename) as src:

        attribs_ = getattr(src, attribs_source_)
        # update meta
        attribs_.update(
            dtype=array_dtype,
            count=1,
            driver="GTiff",
            nodata=nodata_,
            compresion="lzw",
        )

        with rio.open(output_filename, "w", **attribs_) as dst:
            dst.write(array.astype(array_dtype), 1)

check_file = lambda path:  os.remove(path) if os.path.isfile(path) else None

Vector file proc

Incorporamos el archivo vectorial y lo referenciamos en el epsg correspondiente.

gpd_va=gpd.read_file(vector_file)
gpd_va_proj = gpd_va.to_crs(epsg=32720)
#
gpd_vt=gpd.read_file('data/villa_del_totoral.geojson')
gpd_vt_proj = gpd_vt.to_crs(epsg=32720)
#
cases_dict={}

Images proc

Procesamos las imagenes

#
file=zip_20200817
date_='20200817'
dirtmp=tempfile.mkdtemp()
with zipfile.ZipFile(file,"r") as zip_ref:
    zip_ref.extractall(dirtmp)

#
# instanciamos el proc
bp=bands_etl()
# incluimos la capa vectorial
bp.set_mask_from_gpd(gpd_va_proj.geometry)
# Obtenemos los paths
paths,_=bp.get_band_paths(dirtmp, pattern_=['*B0?_10m.jp2','*B1?_20m.jp2'])
#
# Debido a las diferencias en el pixel size de la bandas debemos resamplear
resampled_bands={}
for band in paths:
    case='_'+date_+'_Resampled'
    name=os.path.join(products_tmp,band+case)
    resampled_bands.update({band:name})
    bp.basic_resampling(paths[band],name)
#
shutil.rmtree(dirtmp)
# seteamos las bandas
bp.set_band_names(resampled_bands)
# Generamos los parametros para el recorte
bp.get_clipping_params('B02')
#
case='false_color_urban_'+date_+'.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B12','B11','B04'],tiff_name_=tiff_name_,scale_band=True)
cases_dict.update({case:tiff_name_})
#
case='false_color_urban_'+date_+'_NOTScaled.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B12','B11','B04'],tiff_name_=tiff_name_,scale_band=False)
cases_dict.update({case:tiff_name_})

#
case='rgb_'+date_+'.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B04','B03','B02'],tiff_name_=tiff_name_,scale_band=True)
cases_dict.update({case:tiff_name_})
#
case='rgb_'+date_+'_NOTSCaled.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B04','B03','B02'],tiff_name_=tiff_name_,scale_band=False)
cases_dict.update({case:tiff_name_})

Visualizamos

En falso color y rgb la zona afectada. A modo de referencia incluimos el ejido de villa del totoral (poligono blanco a izquierda).

Ahora vamos a estimar (en forma muy aproximada) el area afectada. Para ello nos valemos de la banda 11. Vemos que es la que mayor contraste parece tener

Analizamos los valores y efectuamos el "thresholding" con el fin de quedarnos con una imagen en blanco y negro.

#collapse
file_path_threshold=os.path.join(products_tmp,'thresholding_b11_20200817.tif')
check_file(file_path_threshold)
!gdal_calc.py --calc="A>1400" -A 'products/false_color_urban_20200817_NOTScaled.tiff' --A_band=2 --outfile={file_path_threshold} --quiet

Vemos que si bien remarcamos la mayoria de la zona, quedan algunas extras. Las intentamos remover e invertimos la imagen.

#collapse
file_path_sieved=os.path.join(products_tmp,'thresholding_b11_sieved_20200817.tif')
check_file(file_path_sieved)
!gdal_sieve.py {file_path_threshold} {file_path_sieved} -st 10

0...10...20...30...40...50...60...70...80...90...100 - done.
<AxesSubplot:>
<AxesSubplot:>

Obtenemos el vector con la zona afectada via gdal-polygonize y visualizamos tanto el raster como el shapefile.

#collapse
vector_fire_sieved=os.path.join(products_tmp,'fire_shape_20200817.shp')
!gdal_polygonize.py {raster_to_shape} {vector_fire_sieved}

0...10...20...30...40...50...60...70...80...90...100 - done.

2020-08-15

Realizamos un procesamiento similar (sin estimación de areas) para una fecha anterior a modo de contraste.

2020-08-20

Observamos la evolución al día 20/08. El procesamiento es similar, en este caso solo extendimos el shape original para enmarcar la zona.

Estimamos nuevamente el area afectada

<AxesSubplot:>

2020-08-22

Observamos la zona afectada al 22/08

#
file=zip_20200822
date_='20200822'
dirtmp=tempfile.mkdtemp()
with zipfile.ZipFile(file,"r") as zip_ref:
    zip_ref.extractall(dirtmp)

#
# instanciamos el proc
bp=bands_etl()
# incluimos la capa vectorial
bp.set_mask_from_gpd(gpd_va_proj.geometry)
# Obtenemos los paths
paths,_=bp.get_band_paths(dirtmp, pattern_=['*B0?_10m.jp2','*B1?_20m.jp2'])
#
# Debido a las diferencias en el pixel size de la bandas debemos resamplear
resampled_bands={}
for band in paths:
    case='_'+date_+'_Resampled'
    name=os.path.join(products_tmp,band+case)
    resampled_bands.update({band:name})
    bp.basic_resampling(paths[band],name)
#
shutil.rmtree(dirtmp)
# seteamos las bandas
bp.set_band_names(resampled_bands)
# Generamos los parametros para el recorte
bp.get_clipping_params('B02')
case='false_color_urban_'+date_+'.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B12','B11','B04'],tiff_name_=tiff_name_,scale_band=True)
cases_dict.update({case:tiff_name_})
#
case='false_color_urban_'+date_+'_NOTScaled.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B12','B11','B04'],tiff_name_=tiff_name_,scale_band=False)
cases_dict.update({case:tiff_name_})

#
case='rgb_'+date_+'.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B04','B03','B02'],tiff_name_=tiff_name_,scale_band=True)
cases_dict.update({case:tiff_name_})
#
case='rgb_'+date_+'_NOTSCaled.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B04','B03','B02'],tiff_name_=tiff_name_,scale_band=False)
cases_dict.update({case:tiff_name_})

Observamos que todavia permanecen algunos focos activos.

Estimamos el area afectada

<AxesSubplot:>

Alternativa para el computo de area afectada

En este caso podemos estimar el area via el NBR (Normalized Burn Ratio). Computamos el indice (_bands_etl().ndi_totiff). Efectuamos el thresholding poligonizamos y generamos el area ( en este caso para el dia 22/08)

(355706.0, 386814.0, 6570194.0, 6620926.0)